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ABSTRACT 

We propose and successfully test against new cosmological simulations a novel analytical 
description of the physical processes associated with the origin of cored dark matter density 
profiles. In the simulations, the potential in the central kiloparsec changes on sub-dynamical 
timescales over the redshift interval 4 > z > 2 as repeated, energetic feedback generates large 
underdense bubbles of expanding gas from centrally-concentrated bursts of star formation. 
The model demonstrates how fluctuations in the central potential irreversibly transfer energy 
into collisionless particles, thus generating a dark matter core. A supply of gas undergoing 
collapse and rapid expansion is therefore the essential ingredient. The framework, based on a 
novel impulsive approximation, breaks with the reliance on adiabatic approximations which 
are inappropriate in the rapidly-changing limit. It shows that both outflows and galactic foun- 
tains can give rise to cusp-flattening, even when only a few per cent of the baryons form stars. 
Dwarf galaxies maintain their core to the present time. The model suggests that constant den- 
sity dark matter cores will be generated in systems of a wide mass range if central starbursts 
or AGN phases are sufficiently frequent and energetic. 



1 INTRODUCTION 

Over the last two decades, galaxies formed in numerical simula- 
tions based on the inflationary ACDM paradigm have suffered from 
a number of well-documented mismatches with observed systems. 
One of the most prominent of these has been the rotation curves 
of disk-dominated dwarf galaxies (e.g. Moore|1994 ; |Flores & Pri- 
mack|1 994; for more recent updates see Simon et aL"l 2005, |Oh et al. 
2011b and references therein). The observed kinematics imply a 
constant density core of dark matter interior to 1 kpc, whereas sim- 
ple physical arguments and simulations suggest that the cold dark 
matter density should be increasing roughly as p « r 1 to vastly 
smaller radii (e.g. |Dubinski & Carlberg|1991 ; Navarro et al. 1996b I. 

Since some of the earliest work on dark matter profiles it has 
been suggested that sufficiently violent baryonic processes might 
be responsible for heating dark matter cusps into cores (Flores 
|& Prim ack 1994). The proposed mechanisms identified in these 
papers fall in two broad categories: supernova-driven flattening 
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Mashchenko et al. 2006 , 2008 ), and dynamical friction from 



infalling baryonic clumps or disk instablities (El-Zant et al.|200"T| 
Weinberg & Katz 2002] |Tonini et al.||2006| [Romano-Diaz et aL] 
20081 |2009| |Pasetto et al.||2010| |Goerdt et al.||2010| |Cole etaij 



201 1[ >. Within the former category, most early works focused on 
a single, explosive mass-loss event. It then became clear that even 
with extreme parameters, such an event transferred insufficient en- 
ergy to dark matter parti cles (|G nedin & Zhao 2002). On the other 
hand, |Read & Gilmore| (2005) showed that several more moder- 
ately violent bursts could be effective in creating a core. Increas- 
ingly sophisticated numerical work by Mashche nko et alT] ( 2006 



|2008| > strongly supported the notion of stellar feedback and energy 
transfer from baryons to dark matter as the generator of cores, but 
did not fully explain the physical mechanism behind this transfer 
or follow the evolution of dwarf galaxies to z = to ensure that the 
cores were long-lived. 

Recently, simulations were able to produce realistic, present- 
day cored dwarf galaxies within a fully cosmological context ( Gov- 
|ernato et al.|2010| henceforth G10). These simulations resolve in- 
dividual star formation 'clumps' at the density of molecular clouds 
leading to galaxies that are additionally realistic because, like 
many observed dwarfs (Dutton 20091, they have no bulge - a con- 
sequence of preferentially expelling low-angular-momentum gas 
from the progenitors via naturally occurring galactic winds ([Brook 
let al.| 201 1; see also Bullo ck et al.| 2001 and jvan den Bosch et al. 
2001). [Oh et al.|p011a[ > confirmed that these effects bring the sim- 
ulations into excellent agreement with observational constraints on 
stellar and Hi content as well as those on overall mass distribu- 
tion. By testing against dark-matter-only runs, G10 provided strong 
support to a model where the core flattening is generated by bary- 
onic effects, in particular by rapid gas motions; moreover, a suite 
of comparison simulations revealed that these effects only become 
significant if stars form in dense clumps (~ 100 amu cm -3 ), sug- 
gesting that energy injection has to be concentrated in local patches 
(see also Ceverino & Klypi nj2~009} . 

The comparison of simulations with recent observations \Oh\ 
|et al.|20 1 la) highlighted that feedback must occur in numerous rel- 
atively mild events to allow thin disks to form. However given the 
lack of analytic framework for understanding the microphysics of 
this process, the precise mechanism of supernova-driven cusp flat- 
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tening was not further elucidated in G10. Providing such a frame- 
work is the aim of the present work. 

The remainder of this paper describes how to model the ef- 
fects of small, central starbursts which create pockets of rapidly 
expanding gas and strong fluctuations in the local potential. Over 
time these repeated processes gradually transfer energy from the 
gas to the dark matter component. This has much in common with 
the view of Mashche nko et al.] ([2006 2008 1 but places more em- 
phasis on disrupting clumps as opposed to pushing them around, 
and clarifies that resonanc^] is not required. Because our picture 
can be modelled mathematically, we are able to validate it against 
the simulations, showing that the envisioned process indeed creates 
cores within the G10 simulations. 

This paper is organized as follows. Section [2] introduces im- 
proved simulations based on those of G10, and discusses the char- 
acteristics of these simulations which predict cusp-flattening, thus 
motivating a study of orbits in rapidly changing potentials (Sec- 
tion [3j. The initial discussion is, for simplicity, limited to power- 
law potentials but Section [4] removes this restriction, presenting 
more general equations to explain the detailed simulation results. 
We relate our work to the wider literature and conclude in Sec- 
tion]^] also discussing the realism of the underlying hydrodynami- 
cal evolution within the simulated galaxies. In a companion paper 
(Governato et al, in prep) we will discuss the scaling of the dark 
matter cores with galaxy masses. 



2 THE SIMULATIONS 

The smoothed particle hydrodynamics (SPH) simulations, run us- 
ing the Gasoline code (Wadsley et al. 2004), are closely related 
to and improve upon those described in more detail by G10. The 
emphasis of the present work is on interpreting dynamical effects, 
rather than discussing the numerical methods in depth; however see 
Section[5]for brief comments on computational accuracy. Our new 
runs output more regular timesteps and include the effects of metal- 
line cooling according to the prescription of Shen et al. (20 10}. The 
simulations in this paper focus on the region hosting the galaxy 
denoted 'DG1' in G10. The 'zoom' technique (e.g. |Katz & White] 
1993 1 allows for a high mass resolution of M p = 3 x IO^Mq (for 
gas particles) and M p = 1 .6 x 1 4 M© (dark matter) with a softening 
of 86 pc in a full ACDM cosmological context. We conducted anal- 
ysis on the two most massive systems within this region: DG1 itself 



3.7 x 10 Mq at z = 0) and a somewhat smaller galaxy 



(M v i r = 1.3 x 1O 1O M at z = 0). Most results will be presented 
for the latter case, because the former undergoes a major merger at 
Z = 3. Although our model does predict the correct flattening for 
DG1, its volatile merger history would introduce undue complexi- 
ties into our discussion. 

In the first run, denoted HT ("high threshold"), stars are al- 
lowed to form only at hydrogen densities exceeding 100 cm -3 . The 
second run, LT ("low threshold"), is identical to the first except that 
it allows stars to form at densities exceeding 0.1cm -3 . 

Adopting the higher threshold for star formation is strongly 
motivated by observational evidence that molecular clouds form at 
such densities (Bigiel et al. 2008; Tacconi et al.||2010) . HT thus 
exhibits more realistic behaviour of the interstellar medium (e.g. 



Saitoh et al. 2008 1. In particular, only at high formation thresholds 
exceeding ~ 10cm -3 will supernova feedback naturally give rise 
to bulk gas motions and outflows (Ceverino & Klypin 20091. This 
follows because individual high-density clumps are efficient in con- 
verting gas to stars, ultimately leading to vast overpressurization of 
the clump from the high local density of supernovae. The particular 
threshold value of 100cm -3 for HT is chosen for consistency with 
G10, in which it was argued that this is the highest density which 
can be considered physical at our present resolution. We have also 
verified that when H2 physics is consistently included in simula- 
tions (Christensen et al, submitted; Governato et al, submitted) star 
formation indeed proceeds only at high densities even in the ab- 
sence of an explicit threshold. 

By contrast, LT's threshold of 0.1cm -3 represents an ap- 
proximate historical norm for galaxy formation simulations (e.g. 
|Navarro & White|1993)|Katz et al.|1996| > broadly compatible with 
the observed cut-off in star formation at low column densities aver- 
aged over ~ kpc scales (e.g. Kennicutt et al. 2007). Until recently 
LT would have been the most motivated choice; however, with the 
addition of metal-line cooling at increasing resolution we now pre- 
fer the HT simulations, using LT as a reference to understand why 
older simulations did not produce the effects of interest here. 

As expected following G10, LT remains cusped, unlike HT 
which develops a 1 kpc dark matter core at z — in both of its two 
most massive halos. In all cases the code consistently follows the 
feedback effects of the stellar populations (Stinson et al. 20061 so 
that, after a delay of ~ lOMyr, significant amounts of thermal en- 
ergy are deposited into the surrounding gas. By z = 2, when HT 
has developed a stable core, LT and HT runs have formed an al- 
most identical mass of stars (7 x 10 7 Mq) and therefore the same 
quantity of supernova energy has been released (7 x lO 56 ergs)0 
The failure of LT to lose its cusp thus reflects a difference in the 
coupling mechanism, not in the absolute energy deposition. 

Figure [T] gives immediate insight into the difference be- 
tween HT (cusp-flattening) and LT (cusp-preserving) simulations 
by showing their spherically-averaged halo density profiles shortly 
before the cusp begins to flatten in HT, at z = 4. Solid, dashed and 
dotted lines indicate respectively dark matter, gas and stellar den- 
sity; thick red and thin blue lines represent the HT and LT runs in 
turn. In the LT run, the gas density cannot exceed the threshold of 
0. 1 cm -3 : stars are able to form and deposit supernova energy, pre- 
venting further cooling. In the HT run, by contrast, the gas density 
rises monotonically towards the centre because most of the gas is 
not eligible to form stars. The result is that the central gas density 
slightly exceeds that of the dark matter. It is natural to suppose from 
this that the halo will have a qualitatively different reaction to the 
presence of baryons in the two runs. (This work focuses on expan- 
sion of dark matter orbits, but we verified that the same processes 
operate on the similarly collisionless star particles in the simula- 
tion, the interesting implications of which are left for future study.) 

Focusing on the HT simulation, Figure [2] (upper panel) shows 
the baryonic mass enclosed within 0.2, 0.5 and 1 .Okpc as a function 
of time. From around 1 .7 Gyr after the big bang, the density near the 
centre undergoes order-of-magnitude fluctuations. First, gas flows 
in, cools and condenses near the centre of the potential well. Then, 
as the density of clumps rises above the 100cm -3 threshold, star 
formation is allowed to proceed. 



1 Although [Mashchenko et"iu~1j2008} described their model as 'resonant' , 
they have since stated that they did not mean to invoke a formal resonance, 
but rather the notion of changes in the potential occurring on roughly the 
dynamical timescale (Wadsley, pri. comm.). 



2 Note, however, that by z = the LT simulation forms 4 x 10 9 M r, in stars 
compared against the HT simulation's 5 x 1() 8 Mq. Only the HT simulation 
forms a realistic dwarf galaxy ^Oh et al.|201 la| . 
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Figure 1. Spherically averaged halo density profiles for high star-formation 
threshold (HT, thick red lines) and low threshold (LT, thin blue lines) simu- 
lations at z = 4, shortly before the dark matter profile starts to flatten in the 
high-threshold model. Solid, dashed and dotted lines show respectively dark 
matter, gas and stellar content. In the low threshold case, dark matter dom- 
inates by orders of magnitude at every radius. In the high threshold case, 
the gas reaches a comparable density to the collisionless matter in the cen- 
tral regions. Gaseous processes can therefore cause heating of collisionless 
components (dark matter and stars) in HT but not LT runs. 



Once supernova energy is dumped into the gas, thermal expan- 
sion forms an underdense bubble of up to several kiloparsecs diam- 
eter (lower panel, Figure[2](. The energy is initially deposited within 
a small volume (a consequence of the high star formation threshold) 
which then reaches temperatures of ~ 10 s K. The gas is vastly over- 
pressurized relative to its surroundings, so expands at close to the 
thermal velocity (typically reaching ~ 300km/s ~ 0.3kpc/Myr). 
Compared to the orbital timescales, which are > 25Myr, the bub- 
ble formation is effectively instantaneous. 

The adiabatic cooling from the expansion is included, but the 
bubbles nonetheless remain hot (~ 10 6 K) by the time they reach 
rough pressure equilibrium with the remainder of the disk. They 
are then sufficiently underdense (~ 10~ 2 cm~ 3 ) that the radiative 
cooling timescale is up to 100 Myr; the bubble can therefore persist 
for this length of time, after which it cools back into the disk if it has 
not actually escaped in a galactic wind. The rapid, repeated fluctua- 
tions in the central mass content of the simulated galaxy are similar 
to those shown in Figure 3 of Mashch enko et al.j ( |2008} . However 
we have verified that, in our case, this is due to the gas being heated 
and expanding outward rather than remaining in rapidly-moving 
coherent clumps as suggested by Mashchenko et al. (2008 1. 

The overall picture for the dark matter is insensitive to the 
ultimate destination of the gas, requiring only the intermittent vari- 
ations explained above; the present work makes no assumptions 
about mass loss. Note, however, that significant winds do exist in 
the simulations; the final baryon fraction in the galaxies is only 
25% of the cosmic value (G10) and only 3% of the baryons have 
been turned into stars. The winds have been shown to be important 
in matching the star formation rates, distribution of stars and final 
baryon fractions of the dwarf galaxies (see McGaugh et al. 2010 
|Brook et al.|20TT]|Oh et al.|201 la| >. 
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Figure 2. (Upper panel) The baryonic mass interior to, from top line to bot- 
tom, 1 kpc, 500 pc and 200 pc (HT simulation). Bursty central star formation 
coupled to strong supernova feedback causes coherent, rapid oscillations in 
the potential interior to 1 kpc. The orbital time of typical dark matter parti- 
cles interior to 1 kpc is > 25 Myr. By contrast the simulated supernova bub- 
bles can encompass the inner kiloparsec in around 3 Myr, far too rapidly for 
the adiabatic approximation to be valid. The lower panel shows the disk- 
plane density during the starburst event at t = 2.56 Gyr, z = 2.67. A large 
underdense bubble has formed at the centre of the disk through thermal ex- 
pansion of gas heated by multiple supernova explosions. The cross marks 
the halo centre. 



3 ANALYTICAL MODEL 

This Section discusses how the energy of a single dark matter par- 
ticle (or star) changes in response to a fluctuating potential sourced 
by gas subject to processes described above. Two restrictions on the 
calculation are imposed throughout the paper: 

(i) the potential is assumed to be spherically symmetric; 

(ii) the tracer particles are assumed to be massless, i.e. the po- 
tential is always external. 

The latter condition represents a decision to focus on the micro- 
physical mechanism via which particles gain energy, rather than 
the subsequent evolution of the self-gravitating system. The first 
condition could in future be relaxed, but makes calculations much 
simpler because particles orbit in the ID effective potential 



V eS (nj,t)=V(r;t) + 



J 

2r 2 ' 



(1) 



where V(r;t) is the time-dependent physical potential and j is a 
conserved angular momentum. 

For simplicity we will temporarily impose two further restric- 
tions which will later be removed in Section[4] 
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(iii) only the normalization of the potential changes, i.e. its func- 
tional form is fixed; 

(iv) the functional form of the potential is a power law. 

Together these imply that the underlying potential in {TJ is specified 
by V (r; f) = Vq (t)r". Some useful test cases fall into this exact form: 
a Keplerian orbit has n = — 1 while a harmonic oscillator implies 
n = 2. The final case will be of particular interest chiefly for its 
analytic simplicity, but we also note that it corresponds to assuming 
a spatially constant density of matter. 

The rate of change of the total energy of a particle orbit- 
ing within the potential, dE /At, is given by the partial derivative 
dV/dt\ r i t y where r(t) denotes the solution to the equations of 
motion. In the limit of an instantaneous change in the potential 
V — > V + AV occurring at time t, the total energy of the particle 
therefore changes by AE = AV (/■(/)) = AVo r{t) n . 

Assuming we have no prior knowledge of the phase of the 
particle, the virial theorem (e.g. Gol dstein et al.|2002| > states that 
the expected potential energy is 



(V} = 



2E 



(2) 



where Eq is the total energy of the particle and the result is inde- 
pendent of /, This and following equations are also therefore valid 
in the one-dimensional (J = 0) subcase. 

If suddenly Vo — ► Vq + AVo, the energy after the potential 
change is Eq + AEi , where 



(A£ 1 }=AV {/'> 



2E AVo 



(3) 



2 + n Vq 

The fiducial adiabatic limit can be attained from here by assuming 
AVo an d AE\ to be infinitesimal and integrating over a series of such 
changes taking Vo smoothly to Vj. This yields a final, finite change 
in energy: 

/VA J /W 

£/,adiabatic = Eq ( — j . (4) 

As expected in the adiabatic limit, equation |4]( implies no energy 
shift, regardless of the intermediate states, if the final potential is 
the same as the initial. This is the central problem of the adiabatic 
approximation. Figure|2]shows that the final distribution of gas will 
indeed be very similar to the initial, and therefore that the adiabatic 
prediction will be for no change in the final distribution of dark 
matter. 

However Section [2] showed that potential changes in the sim- 
ulations take place on timescales much shorter than the dynamical 
time, because the expansion speeds of the supernova-induced bub- 
bles are much larger than the local circular velocity. As the gas 
expands and leaves the galaxy centre, the potential undergoes a se- 
ries of large, instantaneous jumps, invalidating the adiabatic result 
given by equation {4|. 

Instead of integrating, one should therefore recursively apply 
equation {3} to each finite change. If, for instance, the potential 
switches immediately back to its original depth, the second shift in 
energy is given by 

2(£ + <A£ I )) -AV 



(AE 2 ) 



2 + n 



(5) 



Vo + AVo' 

where the angular brackets indicate averaging over the orbital phase 
of our chosen trajectory during both the initial and final instanta- 
neous potential jumps. These conditions are justified because the 
initial blowout is not causally connected to the location of a single 
tracer particle, nor is the exact fractional number of orbits between 



initial blowout and eventual recollapse predictable (this aperiodic- 
ity is illustrated by Figure|2]l. By Taylor expanding, we find that the 
expected final energy of the orbit is given by 



(Ef)=Eo + (AE l +AE 2 }-Eo- 



2n 



Uo j (2+n) 



:E , (6) 



which is always an energy gain for bound orbits since Eq < for 
n < 0. The energy gain is second order in the potential change AVo, 
but linear in the energy. One may verify that, if the potentially first 
changes suddenly but then gradually (i.e. adiabatically) relaxes to 
its original state, we will also see an increase in expected final en- 
ergy of the same magnitude. The essential point is for the initial 
change to be rapid; the energy shift will then follow. 

The special case of the harmonic oscillator (n = 2) is helpful 
in demonstrating the origin of this energy shift, because its dynam- 
ics are especially simple. The potential is separable in Cartesian 
coordinates which means that we can assume the one-dimensional 
subcase without loss of generality. Then, for the case of sudden 
discrete jumps, the analytic form of the solution is written 



x(t) =Acos(cor + t//), 



(7) 



where co 2 (t) = 2Vo(t) while A(t) and yr(t) specify the amplitude 
and phase of the oscillation, which change discontinuously with 
the potential. The new values of A and iff after any jump can be 
determined either through energy arguments as above or by requir- 
ing continuity of both x and x. Without loss of generality, we let CO 
change from Wq to a>i at t = 0. The amplitude of the trajectory after 
the jump, A\, is given by the expression 



cor 



J \ ■ 2 

— sin y/b 



(8) 



which is explicitly dependent on the orbital phase \j/Q at which the 
discontinuity occurs. The corresponding change of energy is 



AE\ = -Eq- 



cot 



- sin y/Q. 



(9) 



We can now analyze the changes in a single trajectory which 
lead us to recover the expected gain in orbital energy under a single 
blowout-recollapse cycle, equation Figure|3](thick line) shows 
an example orbit for which C0q = 10 cof. The initial amplitude is 
unity until, at a certain time, the potential flattens out as mass is lost 
from inside the orbit. Intuitively, or from equation l|9j, this must al- 
ways involve a loss of energy for the particle (since C0q > cof); see 
also panel 1 at the top of Figure [5] However, according to equa- 
tion JS), the new potential is sufficiently flattened that the orbital 
amplitude is now larger than unity (panel 2). This is crucial because 
the energy gain made when the potential change is later reversed (co 
returns from C0i to C0q) will be accordingly larger, scaling with (x 2 ) 
(panel 3). Applying equations {8} and |9j for the reverse jump (i.e. 
with C0\ and C0q exchanged) allow this picture to be directly veri- 
fied. 

This makes more concrete the assertion of equation |6]l that 
one always expects to gain energy during a series of potential 
changes. Yet seen from another perspective, such a claim still needs 
to be reconciled with the underlying dynamics which are fully time- 
reversible. For instance, the time-reverse of Figure |5]( viewing the 
figure from right to left) represents an equally valid trajectory of 
the forwards dynamics and yet loses energy. 

Under time-reversal, the blowout phase maps onto the rec- 
ollapse phase and vice-versa. The irreversibility thus arises not 
from dynamical differences but statistical differences between the 
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Figure 3. The mechanism for injecting energy into the dark matter orbits, 
illustrated by the exact solution for a time-varying harmonic oscillator po- 
tential. The lower panel shows (solid line) a solution to the equations of 
motion where CO 2 = 1 (blue) at early and late times, while at intermediate 
times O) 2 = 0. 1 (red) mimicking baryonic blowout and recondensation. The 
changes in potential occur instantaneously; in this case the final amplitude 
of the oscillation is approximately twice that of the initial orbit. The dashed 
line shows the solution when the potential changes smoothly over several 
orbital periods; this gives adiabatic behaviour, so that the final orbit regains 
the initial amplitude, demonstrating the necessity for relatively sudden po- 
tential jumps. The inset figures (top) illustrate how the post-blowout orbit 
expansion implies that the late-time energy gain dominates over the initial 
energy loss. 



forwards-blowout and reverse-recollapse pictures. In particular a 
uniform prior on the orbital phase before a transition is always as- 
sumed, 2ltp(yn) = 1- On the other hand the phase y/g after the 
transition is determined by 

(On 

tanty/^ = — tanl//Q, (10) 
and, accordingly, the probability distribution function of \/ Q is 



MO 2 J 
<0i 



On 



sin 2 \j/ 



(11) 



The precise functional form l | 1 I [ i is not crucial, only that p(v / o) 
cannot be taken to be uniform. After the sudden baryonic blowout, 
collisionless particles enter their new orbit in a special phase - pref- 
erentially near pericentre - so that they subsequently migrate out- 
wards in unison. 

It is this difference in knowledge of phases before and after 
sudden changes that allows irreversibility in the real universe to ap- 
pear in the model. Only if all collisionless particles were near their 
pericentre just before the baryons returned would the statistical 
properties of the reversed picture match those of the actual model. 
While this is dynamically possible, it is statistically unlikely. 

Finally note that if the changes in potential are introduced 
gradually, the process should become adiabatic and hence re- 
versible. The dashed line in Figure [3] shows a numerical solution 
for which to changes smoothly over several orbital times from too 
to a>i, then back to (On. As expected from equation |4}, the final 
orbital amplitude is the same as its initial value, confirming the 



qualitatively different results to be expected from gradual variation 
as opposed to sudden jumps. 



4 VALIDATING THE ANALYTIC MODEL AGAINST 
SIMULATIONS 

To test the picture expounded above we start by generating a time- 
dependent effective toy potential from the simulations (Section |2j. 
This is given by equation (|TJ, with V(r; t ) calculated from the spher- 
ically averaged density profile. The starting energy £rj and the value 
of j can be determined by specifying initial orbital parameters. The 
angular momentum is necessarily conserved because of the spher- 
ical symmetry of the modelling (restriction 1 of Section[3]l. In the 
simulations the changes in potential are not exactly symmetric (e.g. 
lower panel of Figure [2} ; however we will see below that, for the 
purposes of calculating real-space density profiles, the symmetric 
approximation which enforces constant j actually works extremely 
well. 

As before, the energy shift for one jump is given by averaging 
over possible orbital phases. However the potential V S ph ere is no 
longer an exact power law, so the calculation required is 

J AV eS (r(t))dt 



(AE) 



I 



lit 

AVeff(r)dr 

y/E-Vesir)/ J ^/E-V eH (r) 



dr 



(12) 



where the time integrals are evaluated over an orbital period; af- 
ter changing variables to r this corresponds to integrating over the 
region where the integrand is real. Equation \ \2\ agrees with equa- 
tion {3} for the special case of power-law potentials. 

The remainder of this Section applies expression \\2\ recur- 
sively to a time-series of potentials from the HT (cusp-flattening) 
simulation, at each step updating AV, E and V e ff appropriately^! 
The energy gain is evaluated at every stored simulation timestep; 
the relevant outputs are written every St ~ 27Myr. Thus changes 
occurring on timescales ^ St will implicitly be classified as "rapid" 
(composed of one jump) whereas those occurring on timescales 
3> St will automatically be treated as "adiabatic" (composed of 
many small steps). While the boundary between these limits can- 
not be uniquely defined, the change in behaviour must occur at 
around the orbital period for a particle, which is indeed ~ 25 Myr. 
We verified by running checks with only every second timestep 
(St ~ 54 Myr) that the results presented are insensitive to the pre- 
cise time-slicing. 

The solid lines in Figure[4]show the resulting mean radius (r) 
of orbits as a function of time, where 



(r} = 



/ 



rdr 



yjE{t)-V<s(r,j,t) J ^E(t)~V eS (r;j,t) 



dr 



(13) 



The values of j and En for each orbit are chosen by requiring the 
initial motion to be circular at a range of different radii. As time 
progresses, the orbits starting interior to lkpc migrate outwards, 
reflecting a net gain in energy. Orbits outside this radius are largely 



; Because E takes a random walk, a more accurate result is in principle 
attainable by keeping track of its evolving distribution function rather than 
just its expected value. Our approach here is akin to taking the first term in 
a Fokker-Planck analysis and will be an excellent approximation because 
the energy shifts are approximately linear, as can be shown by generalizing 
equation |5J. 
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Redshift 

4.2 3.9 3.6 3.3 3.0 2.7 2.4 2.1 




1.5 2.0 2.5 3.0 3.5 

t/Gyr 

Figure 4. Using the spherically averaged potential from the simulations, 
we model the expansion of orbits of test particles at different initial radii 
(solid lines). Orbits starting significantly within the inner kiloparsec migrate 
outwards over several gigayears, whereas those starting outside a kiloparsec 
do not feel the rapid potential variations and so remain near their initial 
radius. Our model thus explains the flattening of central density cusps into 
kiloparsec-scale cores in small galaxies through radial outwards migration. 
As expected the reversible, adiabatic model (illustrated for the innermost 
orbit by the dashed line) does not correctly model the heating effect of very 
rapid potential variations in the inner parts of the halo. 



unaffected. In the LT run, by contrast, no tracer particles gain en- 
ergy; those that start on circular orbits, for instance, are predicted 
to remain at the same radius for the entire run. 

Figure [2] implies that the central baryonic potential returns to 
its approximate original shape at the end of each starburst cycle, 
because the gas affected by the supernovae has cooled back into 
the disk (or has flowed out, replaced by fresh gas). In the adiabatic 
limit, where all potential changes occur slowly, the final orbital 
parameters should return to their initial values. Indeed by making 
AV e ff infinitesimal and integrating (12) one obtains 

I y/E(t)- V e ff (r; J, t)dr = constant, (14) 

where again the integral is taken over the real region of the inte- 
grand. This is the generalization of equation Q, and is exactly the 
adiabatic invariant derived through the action-angle approach (e.g. 
|Binney & T remaine 1987). It implies that £g na i = ^initial if the po- 
tential returns to its initial form via a series of slow changes. 

Demanding the adiabatic invariant (14) is constant yields the 
orbital migration in the 'gradual outflows' scenario. The dashed 
line in Figure [4] shows that the result derived in this limit is as ex- 
pected: although temporary changes in the orbital radius do occur, 
they do not persist over time. This underlines the difference be- 
tween our new model (where a tracer particle picks up energy from 
baryons) and the older adiabatic calculations (where the energy of 
a tracer particle is conserved). 

Although Figure [4] shows that orbits gain energy, it cannot be 
used directly to infer the final inner profile of the dark matter. To 
draw conclusions about the evolution of the slope, we evolved the 
energy of ~ 90000 orbits corresponding to all dark matter particles 
in the halo at z = 4. At each timestep, the full radial probability 



.: 2. / :i.r, Gyr 




10" 1 10° 10 1 

r/kpc 



Figure 5. The spherically averaged dark matter density as a function of ra- 
dius, measured at z = 2 when the core has formed in the HT simulations 
(thick dotted line). The solid line shows the density profile at this time 
according to our model (see text for details); this is seen to be in excel- 
lent agreement with the HT simulation. The adiabatic model (dashed line) 
fails to correctly model the cusp flattening, demonstrating the need for the 
improved modelling presented here. The LT comparison simulation (dash- 
dotted line) also remains cusped as explained in Section|2] 

distribution for each particle, 

p(r,E,j)~^==L==, (15) 

was calculated numerically. The sum of the normalized probability 
distributions for all particles then implies a density profile accord- 
ing to 

Pmodei(r) °= -j5Y,p(r;Ei,ji), (16) 

where the sum is over all tracer particles. Time evolution of 
Pmodel ( r ) arises from updating V S phei- e and each at every timestep 
according to equation (12) ; or, for comparison, by solving equa- 
tion (14) to derive the behaviour in the adiabatic limit. 

Starting at z = 4, the distribution function is evolved in this 
way to z = 2; the resulting density profiles are illustrated in Fig- 
ure^ The thick solid line shows our main model [i.e. it is derived 
from equation (12)], and is seen to be in excellent agreement with 
the output of the simulation (dotted line). The dotted line shows 
the results of modelling the baryonic effects using the adiabatic ap- 
proximation [i.e. equation (14)]; the cusp remains, contrary to the 
results of the simulation. This reaffirms that the adiabatic approxi- 
mation does not capture important aspects of the impact of baryons 
on the dark matter. Finally, the dash-dotted line shows the profile 
from the LT (low star formation threshold) simulation which, as ex- 
plained in Section|2] retains its cusp and is therefore in approximate 
agreement with the adiabatically evolved case. 

The calculation described by equation (12) involves calculat- 
ing the particle distribution function for every intermediate step. It 
is possible, therefore, to monitor the rate at which the cusp flat- 
tens and compare it against the simulations. Figure [6] shows the 
time evolution of the measured logarithmic slope at 500 pc for both 
the model density profile, equation (16) , and the simulated density 
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Redshift 
4.2 3.9 3.6 3.3 3.0 2.7 2.4 
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Figure 6. The evolution of the logarithmic slope of the halo in the raw sim- 
ulation data (dotted line) and according to our model (solid line), measured 
at 5()0pc in both cases. Values < — 1 indicate a conventional NFW-style 
cusp; shallower profiles are indicated by values closer to zero. The flatten- 
ing of the profile in the simulations has previously been demonstrated to 
agree with observational data jOh et al.|201 la) . Our physical model for the 
origin of the flattening is in excellent agreement with the detailed simulation 
results. 



profile. As time progresses, both the simulated and model density 
profiles gradually flatten out, with the value rising from < — 1.0 
(cusped) to ~ —0.4, consistent with observations (Oh et al.|2011a| >. 

To conclude, the analytic model presented here predicts flat- 
tening of the slope at the same rate as seen in the HT simulation 
(Figures [5] and [6} ; whereas the adiabatic approximation does not 
predict any significant flattening (Figure|5j. We further verfied that 
the new model did not predict a change in slope for the LT run 
in which the gas density remains too small to generate significant 
fluctuations in the potential. Overall, the new model alone provides 
a convincing explanation for the flattening processes seen by G10. 



5 CONCLUSIONS 

We have proposed a new analytic model that accounts for the flat- 
tening of dark matter cusps into cores. Energy is transferred into 
dark matter particle orbits through repeated, rapid oscillations of 
the central gravitational potential. These oscillations are caused by 
recurrent, concentrated bursts of star formation which induce rapid 
expansion of gas through supernova feedback heating. We veri- 
fied that this process quantitatively accounts for cusp-flattening in 
a novel set of simulations similar to those in G10. The simulations 
include the effects of metal-line cooling, but like those in G10 form 
thin stellar disks and have a galactic star formation efficiency of 
only a few per cent. A comparison simulation (LT) with lower star 
formation density threshold does not form a core, despite forming 
ten times as many stars by z = 0. The model correctly predicts no 
cusp-flattening in this case (Figure|6j confirming our interpretation 
that for cores to form the supernova energy must be injected in a 
concentrated region (Figure[TJ. 

The baryons do not have to escape the system completely, 
but only temporarily vacate the central regions, because the en- 



ergy transfer is inherently irreversible. The G10 simulations do 
exhibit galactic-scale outflows which remove 75% of baryons by 
z = 0. These outflows have other important effects (e.g. Bro ok et al.| 
|20TT] >, yet the mass involved is only around a third of the mass 
involved in the central blowout-recollapse cycle. The relation be- 
tween the galactic outflows and the local feedback will be the focus 
of future work. 

The picture is related to the more established view that re- 
moval of baryons through galactic superwinds could cause the cen- 
tral dark matte r profile to flatten (e.g. | Navarro et al.|1996a||Gnedin| 
|& Zhao|200"2l |Read & Gifmore|2005) . However, our model con- 
firms that extreme, violent mass-loss events are not necessary, as 
suggested by (Mashchenko et al. 2006, 2008 I. This is important be- 
cause the more moderate heating events allow retention of baryons 
and the formation of a thin stellar disk. 

Dynamical friction from infalling baryonic clumps (El-Zant 



eTaIl[20071 [Mo~& Mao 2004) |Romano-Di'az et al.|[2009l TGoerdt 



et al.|20 10l does not appear to play a dominant role in our simu- 



lations. In the LT simulations, no cores form; therefore effects of 
dynamical friction are ruled out except from the densest clumps in 
HT. In the HT simulations, any dense infalling clumps are typically 
disrupted long before they reach the inner kiloparsec. To test the 
effect of dynamical friction, one would first need to remove the ex- 
plosive events associated with feedback. However as the feedback 
energy is decreased, dense clumps start to pile up in galaxies until 
time integration becomes computationally unfeasible, at the same 
time creating an unrealistic dense central bulge. Additionally we 
do not have the resolution to produce star clusters iGoer dt et al.| 
|2010| > which could be more robust to disruption than gas clumps. 
As a result we are not ruling out dynamical friction as an agent for 
weakening cusps except for the particular simulations in use here. 

The DG1 case (which undergoes major mergers at z = 3 and 
z = 1; see Section [5} has a very similar cusp-flattening history to 
the galaxy described in this work. The mergers themselves have 
no measurable effect on the halo slope and, as expected through 
arguments based on Liouville's theorem, cores remain present after 
the merger (Kazantzidis et al. 2006). In Governato et al (in prep.) 
we show that cores form in the large majority of dwarf galaxies, 
irrespective of their assembly history. 

The model is reminiscent of the violent relaxation envisioned 
by |Lynden^ Bell ( 1967), although our analysis does not attempt to 
generate the gravitational potential self-consistently. It would thus 
be desirable in the future to work the microphysics into a broader, 
analytical description of the evolution of a self-consistent dark mat- 
ter distribution function. We would further like to investigate the 
importance of departures from exact spherical symmetry. The sim- 
ulated supernova explosions are rarely exactly on the axis of the 
disk, so even axisymmetry is violated. This will lead to the non- 
conservation of angular momentum, which must be understood be- 
fore we can investigate the impact of the process on the anisotropy 
of the orbits (Tonini et al. 2006). 

Stars, like dark matter particles, are collisionless and therefore 
should be subject to the same migratory processes outlined here. 
Stars forming near the centre of the HT galaxies indeed migrate 
outwards; it could be natural in this context that the scale-length of 
the dark matter cores and the stellar disks are approximately equal, 
an observational relation noted by Ge ntile et al.| p009). Indeed the 
scale-lengths of both stellar disk and dark matter core of the dwarf 
galaxies discussed here are approximately 1 kpc (G10; Broo ks et al.| 
201 1). To make this link convincing will require a more systematic 
study of scaling with mass (Brooks et al in prep). 

While the analytic model we have described is independent of 
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the detailed gas dynamics, to obtain results the simulated hydro- 
dynamics are used as an input. The cusp-flattening effect is there- 
fore only achieved in reality if the rapid gas motions predicted by 
the SPH code are reproduced. Mesh-based codes (Teyssier 2002 
|Q'Shea et aL]|2004| [Springel 2010) highlight potential shortcom- 
ings of traditional SPH such as its poor handling of instabilities 
related to sharp density contrasts ([Agertz et al. 2007] |Bauer &] 
Springel 20111. In future we intend to address the sensitivity of 
our results to these inaccuracies through comparison with alterna- 
tive codes and use of forthcoming improvements to the Gasoline 



SPH engine that reduce artificial surface tension (see also |Read & 
Hay field 201 1 1. Through direct comparison of various codes Scan- 
napieco et al. |20TT| conclude that, for most practical purposes 



the choice of sub-grid model approximation (via which supernova 
energy is coupled to the gas) is more critical than the numerical 
technique. It will be of interest to determine whether other feed- 
back mechanisms and numerical methods can reproduce our re- 
sults. Indeed recently Martizzi et al. (2011) reported the formation 
of dark matter cores in AMR simulations through gas fluctuations 
very similar to ours, but driven by AGN activity in clusters. 

We are currently investigating how the cores scale for 
1O 9 M < M vir < 1O 12 M (Governato et al, in prep). At higher 
masses the results will depend on a detailed interplay between the 
deepening dark matter potential, increased star formation rates, and 
the nature of AGN feedback (Martizzi et al.|20lT} . The challenge 
of running suitable simulations to tackle the most massive systems 
at sufficient resolution is formidable, but one that we hope to tackle 
in due course. 
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